Consignes

Complétez ce document en remplissant les chunks vides pour écrire le code qui vous a permis de répondre à la question. Les réponses attendant un résultat chiffré ou une explication devront être insérés entre le balises html code. Par exemple pour répondre à la question suivante :

La bioinfo c'est : <code>MERVEILLEUX</code>.

N’hésitez pas à commenter votre code, enrichier le rapport en y insérant des résultats ou des graphiques pour expliquer votre démarche. N’oubliez pas les bonnes pratiques pour une recherche reproductible ! Nous souhaitons à minima que l’analyse soit reproductible sur le cluster de l’IFB.

Les informations essentielles que nous attendions étaient la version des outils utilisés et le détail des commandes lancées. L’organisation des répertoires est propre à chacun, l’important étant d’organiser les fichiers de façon intelligible.

Introduction

Vous allez travailler sur des données de reséquençage d’un génome bactérien : Bacillus subtilis. Les données sont issues de cet article :

Analyses

Organisation de votre espace de travail

mkdir ~/EVALUATION
cd ~/EVALUATION
mkdir RAW_DATA REFERENCES QC CLEANING MAPPING

Téléchargement des données brutes

  • Récupérez les fichiers FASTQ issus du run SRR10390685 grâce à l’outil sra-tools [1]
module load sra-tools/2.10.3
fasterq-dump --version
# "fasterq-dump" version 2.10.3
srun --cpus-per-task 8 fasterq-dump -S -p SRR10390685 --outdir RAW_DATA --threads 8
cd RAW_DATA
gzip SRR10390685_1.fastq
gzip SRR10390685_2.fastq 
cd -
  • Combien de reads sont présents dans les fichiers R1 et R2 ?
cd ~/EVALUATION/RAW_DATA

module load seqkit/0.14.0
seqkit stat *.fastq.gz > seqkit.tsv
head seqkit.tsv

# file                    format  type   num_seqs        sum_len  min_len  avg_len  max_len
# SRR10390685_1.fastq.gz  FASTQ   DNA   7,066,055  1,056,334,498       35    149.5      151
# SRR10390685_2.fastq.gz  FASTQ   DNA   7,066,055  1,062,807,718      130    150.4      151

# other solution but less generic

# zcat SRR10390685_1.fastq.gz | echo $((`wc -l`/4))
# zcat SRR10390685_2.fastq.gz | echo $((`wc -l`/4))

Les fichiers FASTQ contiennent 7 066 055 reads.

  • Téléchargez le génome de référence de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse
cd ~/EVALUATION/REFERENCES
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.fna.gz
gzip -d GCF_000009045.1_ASM904v1_genomic.fna.gz
  • Quelle est la taille de ce génome ?
cd ~/EVALUATION/REFERENCES
module load seqkit/0.14.0
seqkit version
# seqkit v0.14.0

seqkit stat GCF_000009045.1_ASM904v1_genomic.fna
# file                                  format  type  num_seqs    sum_len    min_len    avg_len    max_len
# GCF_000009045.1_ASM904v1_genomic.fna  FASTA   DNA          1  4,215,606  4,215,606  4,215,606  4,215,606

La taille de ce génome est de 4 215 606 paires de bases.

  • Téléchargez l’annotation de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse
cd ~/EVALUATION/REFERENCES
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.gff.gz
  • Combien de gènes sont connus pour ce génome ?
zgrep -v "^#" GCF_000009045.1_ASM904v1_genomic.gff.gz | awk '($3 == "gene")' |wc -l
# 4448

4 448 gènes sont recensés dans le fichier d’annotation.

Contrôle qualité

  • Lancez l’outil fastqc [2] dédié à l’analyse de la qualité des bases issues d’un séquençage haut-débit
cd ~/EVALUATION/QC
module load fastqc/0.11.9
fastqc --version
# FastQC v0.11.9

for i in ../RAW_DATA/*.fastq.gz ; do srun --cpus-per-task 8 fastqc $i -o . -t 8 ; done

module load multiqc/1.9
multiqc --version
# multiqc, version 1.9

srun multiqc -d . -o .
  • La qualité des bases vous paraît-elle satisfaisante ? Pourquoi ?

La qualité des bases me paraît satisfaisante car la qualité moyenne des bases est supérieure à 30 sur toute la longgueur des reads comme le montre le graphique Sequence Quality Histograms ou le graphique Per Sequence Quality Scores. De plus, il y a très peu de bases indéterminées (N).

Lien vers le rapport MulitQC

  • Est-ce que les reads déposés ont subi une étape de nettoyage avant d’être déposés ? Pourquoi ?

Oui, car la distribution de taille des reads (Sequence Length Distribution) montre que certains reads sont plus petits que 151, la longueur attendue. Néanmoins, cela représente très peu de reads donc ce n’est pas vraiment problématique. Il serait intéressant de trouver l’information des filtres appliqués (via la publication pour des données publiques ou via la plateforme de séquençage pour des données non publiées.

  • Quelle est la profondeur de séquençage (calculée par rapport à la taille du génome de référence) ?
Formule à appliquer :

(length all R1 + length all R2) / Taille du génome

soit ici:

(1056334498+1062807718)/4215606

La profondeur de séquençage est de : 502.6898187 X.

Nettoyage des reads

Vous voulez maintenant nettoyer un peu vos lectures. Choisissez les paramètres de fastp qui vous semblent adéquats et justifiez-les.

cd ~/EVALUATION/CLEANING

module load fastp/0.20.0
fastp --version
# fastp 0.20.0

srun --cpus-per-task 8 fastp --in1 ../RAW_DATA/SRR10390685_1.fastq.gz --in2 ../RAW_DATA/SRR10390685_2.fastq.gz --out1 SRR10390685_1.fastq.gz --out2 SRR10390685_2.fastq.gz --html fastp.html --thread 8 --cut_mean_quality 30 --cut_window_size 8 --length_required 100 --cut_tail --json fastp.json

seqkit stat SRR10390685_[12].fastq.gz > seqkit.txt
head seqkit.txt

# file                    format  type   num_seqs      sum_len  min_len  avg_len  max_len
# SRR10390685_1.fastq.gz  FASTQ   DNA   6,777,048  996,891,051      100    147.1      151
# SRR10390685_2.fastq.gz  FASTQ   DNA   6,777,048  990,442,597      100    146.1      151

Les paramètres suivants ont été choisis :

Parametre Valeur Explication
–cut_mean_quality 30 pour un score moyen dans la fenêtre glissante > 30
–cut_window_size 8 pour une taille de fenêtre glissante de 8
–length_required 100 pour ne garder que les reads de taille > 100
–cut_tail pour faire partir la fenêtre de l’extrémité 3’ du read

Ces paramètres ont permis de conserver 6 777 048 reads pairés, soit une perte de 4.0900757 % des reads bruts.

Alignement des reads sur le génome de référence

Maintenant, vous allez aligner ces reads nettoyés sur le génome de référence à l’aide de bwa [3] et samtools [4].

cd ~/EVALUATION/MAPPING

module load samtools/1.10
samtools --version
# samtools 1.10
# Using htslib 1.10.2

module load bwa/0.7.17
bwa
# Version: 0.7.17-r1188

srun bwa index ../REFERENCES/GCF_000009045.1_ASM904v1_genomic.fna

srun --cpus-per-task=4 bwa mem ../REFERENCES/GCF_000009045.1_ASM904v1_genomic.fna ../CLEANING/SRR10390685_1.fastq.gz ../CLEANING/SRR10390685_2.fastq.gz -t 3 | samtools view -hbS - > SRR10390685.bam
srun samtools flagstat SRR10390685.bam > SRR10390685.bam.flagstat
srun samtools sort SRR10390685.bam -o SRR10390685_sorted.bam
srun samtools index SRR10390685_sorted.bam
  • Combien de reads ne sont pas mappés ?
samtools view -f 4 -c SRR10390685.bam
# 744540

744 540 reads ne sont pas mappés.

Croisement de données

Calculez le nombre de reads qui chevauchent avec au moins 50% de leur longueur le gène trmNF grâce à l’outil bedtools [5]:

Pour répondre à cette question, je choisis de récupérer tous les reads qui chevauchent le gène en faisant l’intersection des reads du fichier BAM avec le fichier GFF3. Ensuite je compte le nombre de reads du fichier BAM généré.

cd ~/EVALUATION/MAPPING

module load bedtools/2.29.2
bedtools --version
# bedtools v2.29.2

zgrep trmNF ../REFERENCES/GCF_000009045.1_ASM904v1_genomic.gff.gz | awk '$3=="gene"' > trmNF.gff3
bedtools intersect -a SRR10390685_sorted.bam -b trmNF.gff3 -f 0.5 > SRR10390685_on_trmNF.bam
samtools view -c SRR10390685_on_trmNF.bam
# 2801

2 801 reads chevauchent le gène d’intérêt.

Visualisation

Utilisez IGV [6] sous sa version en ligne pour visualiser les alignements sur le gène. Faites une capture d’écran du gène entier.

Pour répondre à cette question, il faut récupérer en local le fichier BAM et son index (.bam.bai) ainsi que le génome de référence, et éventuellement le GFF3 du gène d’intérêt pour observer les parties flanquantes.

IGV : couverture du gène trmNF

Synthèse

tree ~/EVALUATION
EVALUATION
├── CLEANING
│   ├── fastp.html
│   ├── fastp.json
│   ├── seqkit.txt
│   ├── SRR10390685_1.fastq.gz
│   └── SRR10390685_2.fastq.gz
├── MAPPING
│   ├── SRR10390685.bam
│   ├── SRR10390685.bam.flagstat
│   ├── SRR10390685_on_trmNF.bam
│   ├── SRR10390685_sorted.bam
│   ├── SRR10390685_sorted.bam.bai
│   └── trmNF.gff3
├── QC
│   ├── multiqc_data
│   │   ├── multiqc_data.json
│   │   ├── multiqc_fastqc.txt
│   │   ├── multiqc_general_stats.txt
│   │   ├── multiqc.log
│   │   └── multiqc_sources.txt
│   ├── multiqc_report.html
│   ├── SRR10390685_1_fastqc.html
│   ├── SRR10390685_1_fastqc.zip
│   ├── SRR10390685_2_fastqc.html
│   └── SRR10390685_2_fastqc.zip
├── RAW_DATA
│   ├── seqkit.tsv
│   ├── SRR10390685_1.fastq.gz
│   └── SRR10390685_2.fastq.gz
└── REFERENCES
    ├── GCF_000009045.1_ASM904v1_genomic.fna
    ├── GCF_000009045.1_ASM904v1_genomic.fna.amb
    ├── GCF_000009045.1_ASM904v1_genomic.fna.ann
    ├── GCF_000009045.1_ASM904v1_genomic.fna.bwt
    ├── GCF_000009045.1_ASM904v1_genomic.fna.pac
    ├── GCF_000009045.1_ASM904v1_genomic.fna.sa
    └── GCF_000009045.1_ASM904v1_genomic.gff.gz

6 directories, 31 files

References

1. toolkit NS. NCBI sra toolkit. NCBI, GitHub repository. 2019.

2. Andrews S. FastQC a quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.

3. Li H. Aligning sequence reads, clone sequences and assembly contigs with bwa-mem. arXiv preprint arXiv:13033997. 2013.

4. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and samtools. Bioinformatics. 2009;25:2078–9.

5. Quinlan AR, Hall IM. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.

6. Thorvaldsdóttir H, Robinson JT, Mesirov JP. Integrative genomics viewer (igv): High-performance genomics data visualization and exploration. Briefings in bioinformatics. 2013;14:178–92.

 

A work by Migale Bioinformatics Facility

https://migale.inrae.fr

Our two affiliations to cite us:

Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France

Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France